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Abstract. We formulate and solve a model of dynamical arrest in colloids. A particle is coupled to the bath 
of statistically identical particles. The dynamics is described by Langevin equation with stochastic external 
force described by telegraphic noise. The interaction with the bath is taken into account self-consistently 
through the back-reaction mechanism. Dynamically induced glass transition occurs for certain value of 
the coupling strength. Edwards- Anderson parameter jumps discontinuously at the transition. Another 
order parameter can be also defined, which vanishes continuously with exponent 1/2 at the critical point. 
Non-linear response to harmonic perturbation is found. 
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1 Introduction 

Glass transition and slow relaxation in systems charac- 
terised by weak ergodicity breaking remains still an open 
area, despite many efforts and numerous significant re- 
sults in the last decade [T|P|nj||I||^||^|[7||Bl Plinillllll2lllHI 

phenomena we are motivated here mostly by the effect 
of dynamical arrest in colloidal matter [2311241125112611271 
I28II29II30II31| , observed experimentally and thoroughly in- 
vestigated by numerical simulations and Mode-Coupling 
method. Below the transition point, the dynamics effec- 
tively leads to a glassy state with diverging viscosity, how- 
ever the static thermodynamic transition may not be iden- 
tifiable. Indeed, dynamical or structural arrest demon- 
strates the glass transition as a purely dynamic and self- 
consistent phenomenon, where casual slowdown of cer- 
tain particles prevents some other particles from moving, 
which may slow down the others even more etc. The self- 
consistent nature of the phenomenon is reflected by the an- 
alytical approaches available now. One of the most striking 
phenomena in colloids, suspensions and granular matter is 
the non-Newtonian response to mechanical perturbation. 
On one hand, we can have shear-thinning, which amounts 
to a decrease if viscosity due to applied field, which can 
be interpreted as restoration of ergodicity due to pertur- 
bation |25| . On the other hand, increase of viscosity may 
result in shear thickening or even jamming, typically ob- 
served in particulate or granular matter |32II33| . 
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From the theoretical side, the Mode-Coupling (MC) 
equations provide us with a well-established framework, 
capable of explaining a good deal of experimental data 
|Mll35ll36llTf] . The attempts to derive the MC equations 
starting from the Hamiltonian of the system were suc- 
cessful in the mean-field approximation. It was perhaps 
the p-spin spherical model P1I51Id1|13| . where the machin- 
ery reached the farther edge of our current understanding 
of the phenomenon. 

However, the bottom-up approach starting with writ- 
ing explicit Hamiltonians is far from being complete. The 
presence of the reparametrisation invariance [131138) leaves 
the numerical solution of the MC equations as the only 
means for obtaining the true time-dependence of the cor- 
relation and response functions. Also the mean-field ap- 
proximation generally used now seems to be very difficult 
to overcome. 

The serious difficulties remaining in using the more 
advanced MC techniques leave the space for more simple 
phenomenological approaches. We want to follow this path 
in the present work. 

Indeed, the mathematical substance of the Mode Cou- 
pling method can be summarised by saying that the time 
dependence of the correlation (and response) functions de- 
pends non-linearly and in time-delayed manner on these 
functions themselves. Actually, the memory kernel in the 
MC equations, which is primarily dictated by the proper- 
ties of the reservoir, depends of the system dynamics. 

We may represent the dynamics of the system by a 
stochastic process and the parameters of the process de- 
pend on time through the averaged properties of the pro- 
cess itself. In order to study generic properties of such 
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problems it can be useful to establish a simple idealised 
model, which would capture the essential mathematic in- 
gredients while avoiding the complications which arise from 
choosing a specific Hamiltonian at the beginning. The 
most important ingredient in such an idealised model should 
be the mechanism of the back-reaction. 

We introduced recently a very simple stochastic 
process, in which the back-reaction leads to rich dynamic 
behaviour. The main characteristics was the presence of a 
phase transition from ergodic to non-ergodic phase. The 
principal aim of the present work is to investigate analyti- 
cally some of the properties of the transition and from the 
numerical solution of the corresponding differential equa- 
tions infer the non-trivial critical behaviour. 



2 Langevin equation with back-reaction 
2.1 System of coupled particles 

The model system we will have in mind will be composed 
of particles, relaxing to their equilibrium positions under 
the influence of surrounding particles. They can be viewed 
as colloidal particles immersed in a solvent, but the for- 
mulation of our model is generic enough to allow for other 
interpretations as well, e. g. they can be viewed as micro- 
domains in a relaxor ferroelectric material. 

The time evolution of the model can lead to dynami- 
cal arrest, where particles are locked in their positions by 
surrounding particles, which are also locked in their turn. 
Therefore, the dynamics can lead to the spatially disor- 
dered but time-stable stationary state with glass proper- 
ties. The indication of the glass transition will be the non- 
zero value of the Edwards- Anderson order parameter and 
sensitivity to initial conditions. The interaction between 
particles will be taken into account on a phenomenologi- 
cal level; if we concentrate on a randomly selected parti- 
cle (single relaxor), the external field from the rest of the 
system (reservoir) will change as the states of the other 
particles (relaxors in reservoir) evolve. The changes in the 
local external field will be the more rapid the faster is the 
evolution of the other particles. This leads to the idea of 
expressing the intensity of the changes in th external field 
through the velocity of movement of the relaxors in the 
reservoir. As we suppose all particles to be statistically 
identical, the movement of our single relaxor should be 
in probabilistic sense equivalent to the movement of any 
relaxor within the reservoir. This consideration closes the 
loop. 

We will try to express the intensity of the changes 
of the local field through the averaged properties of the 
movement of the single relaxor itself. This introduces the 
idea of back- reaction: the probabilistic properties of the 
reservoir dictate the system evolution and the averaged 
system dynamics tunes the properties of the reservoir it- 
self. 

To be more specific, our single relaxor will be described 
by the continuous real stochastic variable X{t). It will 
evolve under influence of the environmental force, repre- 
sented by the stochastic variable Q{t). The force will be 



modelled by a two-valued random process Q{t) E {—q, +q} 
|4()j . jumping at random instants. Occurrence of the jumps 
are governed by a self-exciting point process |41II42| with 
time-dependent intensity ^X{t). 

For given (friction-reduced) force Q{t) the single re- 
laxor is described by the Langevin equation [^44II45| 

|x(t) = -7X(i) + Q(t) (1) 

with initial conditions X(0) = Xq and Q(0) = Qq- 

We may consider the process X(t) as a movement of an 
over-damped particle which slides in a parabolic potential 
well; the parabola jumps between two positions at random 
instants, with time-dependent rate ^A(i). 

The intensity of the process, or the frequency of jumps, 
^A(i), is related to the movement of the surrounding re- 
laxors, considered as a reservoir. Intuitively the frequency 
must be smaller if the movement of the relaxors is slower. 
Therefore, the function X{t) should be coupled to the ve- 
locity V(t)= ^^ X(t). 

However, it is not obvious a priori what should be the 
specific functional dependence. We only require that the 
dependence is described by a non-negative function an- 
alytic at the origin. The simplest choice satisfying this 
property is 

A(t) = 4(V2(i)) (2) 

where e is the dimensionlcss coupling strength parameter. 
The latter prescription is the form of the back-reaction we 
will study in the following. The model described above is 
sketched schematically in the Fig. ^ 




£x(0=-YX(r) + Q(0 

Fig. 1. Schematical picture of our model. Relaxing particles in 
the reservoir influence the relaxation of the selected particle, 
described by the Langevin equation shown in the frame box. 

The parameters 7 and q can be in principle rescaled 
to 1 by appropriate choice of the units of time and length. 
Therefore, the coupling strength e remains to be the only 
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physically relevant parameter tuning the behaviour of the 
system. As we will see, there is a qualitative change in the 
behaviour of the system at a certain critical value of e. 



2.2 Properties of the environmental force Q(t) 

The force Q{t) is a time-inhomogeneous Markov process, 
therefore its properties are fully described by a master 
equation. More specifically, let us define the probabilities 

7r±(t) =Prob{Q(t) =±<z} (3) 



making the vector 7r(t) = 




which satisfies the 



Pauli master equation 

±„it)^-^\{t)(^_\-l^nit). (4) 

Solving the equation amounts to calculation of the corre- 
sponding time-ordered exponential. The averages and cor- 
relation functions can be expressed through the integrated 
intensity 

A{t) = f X{t')dt' . (5) 
Jo 

The time dependence of the vector 7r(i) can be written 
through the semi-group operator 

««.«„)=(j?)+i(e-<'>«""'-i)(_;-;) (6) 

as 7r(t) = R{t, to)TT{to). Note that the semi-group operator 
obeys R{t,T)R{T,to) = R{t,to) Vr S [tQ,t] which testifies 
the Markov property of the process. 
More explicitly, we find 

(Qit)) = (Qo) cxp(-^(i)) (7) 



{Q{t)Q{h)}^q'cxp{-\A{t)~A{t,)\) (8) 

Similarly, also the higher correlation functions can be writ- 
ten as products of exponentials with combinations of A{t) 
with appropriate time arguments in the exponents. In fact, 
higher order correlation functions factorise into product of 
first and second order correlations, e. g. 

2k k 

{l[Q{U))=ll{Q{hi)Q{hi-i)) (9) 

1=1 1=1 

for t2k > t2k-i > ••• > ^2 > ii- Another consequence is, 
that the cumulants of higher order than two vanish. 

Note that the function A{t) is non-decreasing and can 
either diverge (if limj^cxj A(t) > 0) or assume a finite limit 
for t oo, if X{t) approaches fast enough. 
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3 Glass transition and asymptotic relaxation 
3.1 Equations for moments 

For any given realisation of the process Q{t), the formal 
solution of Eq. |^ is 

X(t) =Xoc-T*+ r e-T(*-*')Q(t')dt' . (10) 

If the function X(t) were known, various moments (and 
correlation functions) of the random process X(t) could 
have been computed from H1U|) using the expressions Q 
and However, in our case the function X{t) should be 
computed from the condition (jSJ, relating it to the second 
moment of the time derivative of X(i). This suggests that 
sufficiently broad set of moments of X{t) and V(t) may 
provide a closed set of ordinary differential equations. The 
solution of this set will yield the closed description of the 
behaviour of our model. 

Indeed, we can define four auxiliary functions 

si(t) = e-^(*) (11) 

S2{t) = e^-^* / dt'c^*'-^'^*'^ (12) 
Jo 

S3(i) = e-^*-^(*) /*dt'e^*'+^(*') (13) 
Jo 

S4(t) = e-2'r* rdi'e^*'-^(*') /* di" e^*"+'^(*") (14) 
Jo Jo 

and express the requested quantities through these func- 
tions. For example the average coordinate can be written 
as 

(X(t)) = (Xo)c-^* + (Qo)s2(t). (15) 
Similarly, the second moment of the coordinate is 

(X2(<)) = {XD e-27t+2(XoQo) e"^* ■S2{t)+2q^ s^it) . (16) 

The functions si{t) to S4{t) can be found by solving the 
set of non-linear differential equations 



h{t) = -X{t)sj{t) (17) 

S2it) ^ -1 S2it) + Si{t) (18) 

i3(t) = l-(7 + A(t)) S3(t) (19) 

hit) ^ s^it) + ssit) (20) 



with initial conditions si(0) = 1, S2(0) — 33(0) ~ 54(0) = 
0. The function A(i) occurring in the latter equations is 
itself a combination of the functions si{t) to 84(1) 

A(t) =e[7-272 S3(i)+273.s4(i)] + 

+ '-f biK) e"'^* + 27(XoQo) e--^* 82(1)- (21) 

-2(XoQo)e^''*si(t)] . 

In the following we will assume that the initial con- 
dition of the stochastic process is Xq = and (Qo) = 0, 
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Fig. 2. Time evolution of the auxiliary functions si{t) (red 
line), S2{t) (blue line), ssit) (green line), and S4(t) (magenta 
line) for 7 — 1, q = 1. The panel a) corresponds to the value of 
the parameter e = 1.2, while in the panel b) we have e — 0.8. 



Fig. 3. Time evolution of the average coordinate (X(t)) com- 
puted for initial condition (Qo) — 1 (red line) and of the the 
switching rate X{t) (blue line) for 7 = 1, g = 1. The panel a) 
corresponds to the value of the parameter e — 1.2, while in the 
panel b) we have e = 0.8. 



except explicitly mentioned cases. This lead to significant 
simplification of the mathematical structure of the equa- 
tions. Indeed, the equations for S3{t) and S4(t) form a 
closed pair of equations 

S3(i)=l- (1 + 6)753(0 + 

+ 2e72(s3(i)-7S4(t))s3(i) (22) 
hit) =S3{t) ~ 2j Siit) . 

Unfortunately, the system H22(l cannot be solved analyti- 
cally. The best one can do is to transform the non-linear 
Ricatti-type set ()22|l to one differential equation of Abel 
type, whose solution, however, is not generally known. 
Therefore, we will solve the equations (|22|l numerically. 
Nevertheless, there is still a significant amount of infor- 
mation which can be extracted analytically. 

The typical results of numerical solution are shown in 
figures 13 and In Fig. |5] we can see the time evolution 
of the auxiliary functions si{t) to S4{t). Fig. 13 shows the 
evolution of the switching rate A(t) and average coordinate 
(X(t)) for non-zero value of the initial condition (Qo). We 



can observe qualitatively different behaviour for e < 1 
and e > 1: first, the switching rate approaches a non- 
zero limit for e > 1, while for e < 1 it decays to zero. 
This means that in the latter case the system effectively 
freezes. This is further confirmed by the observation that 
for e < 1 the limit value of the average coordinate depends 
on the initial conditions, while in the opposite case the 
dependence on initial conditions is lost for large times, the 
system equilibrates and the average coordinate converges 
always to zero. The following sections are mainly devoted 
to the analytical investigation of the above observations. 

3.2 Fixed points 

The first step in investigating the behaviour of the sys- 
tem H22|l is the search for the fixed points [sg , sj] of the 
dynamics. We found that there are only two fixed points, 
namely 

14,4] ^[-,-^] (23) 
76 27^6 
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and 



(24) 



Let us denote Aoo the value of X{t) calculated at the corre- 
sponding fixed point. Using H2f|) the fixed point 123f) yields 
Aoo = (e — 1)7 and the fixed point yields Xoo = 0. 

The linear stability analysis reveals that for e > f the 
fixed point H23|) is stable, while (|24|l is unstable. On the 
other hand, for e < f the fixed point H24(l is stable, while 
H23(l is unstable. The case e = f is a marginal one, where 
both fixed points have one of the eigenvalues equal to 0. 
Therefore, the value e = f marks a transition, whose na- 
ture will be further pursued in the following. 



where x = x^/q, 0{a) is the Heaviside unit-step function, 
and B{a, b) denotes the Beta-function We can observe 
a qualitative change at the value e = 3. For e > 3 the 
limiting distribution H28(l has a maximum for x ~ and 
approaches at the edges of the support [— g/7, q/j], while 
for e < 3 it has a minimum at x = and diverges at the 
edges of the support. The tendency for accumulating the 
probability close to the points ±(7/7 when e decreases can 
be regarded as a precursory phenomenon of the transition 
to the non-ergodic regime, investigated in the next sub- 
section. 



3.4 Non-ergodic regime e < f 



3.3 Ergodic regime e > 1 

In this case the relevant fixed point is (|23() and inserting 
its value to the expressions for the moments of X(t) we find 
that both the average coordinate and the average velocity 
relaxes to zero. On the other hand, the fluctuations of the 
coordinate reach positive value, so 



lim {X{t)) = 

t—*oo 



hm (X^W)-^ 

i— i-oo 'ye 



(25) 



The external force switching rate converges to positive 
constant lim4_,oo X(t) = 7 (e— I). In all cases the quantities 
of interest converge exponentially to their limit values The 
rate of convergence is determined by the lowest in absolute 
value eigenvalue, which is 



Ml 



(26) 



In the interval e G (4 — 2 V2, 4 + 2\/2) the eigenvalue ac- 
quires a non-zero imaginary part, which means that os- 
cillatory behaviour is superimposed over the exponential 
relaxation. 

The overall picture is the following. The back-reaction 
leads to self-adjustment of the switching rate of the ex- 
ternal force exerted by the reservoir. The coordinate fluc- 
tuates around the origin and these fluctuations are sta- 
tionary. Therefore, the stationary regime of the system 
corresponds to the primitive version with fixed A, except 
the fact that the value of A is not given from outside, but 
tuned by the dynamics itself. We call this regime ergodic, 
because the particles do not freeze at some value of the 
coordinate X(t) but fluctuate forever. 

The probability density for the coordinate 



P(x,t) = — Prob{X(t) < x} 
ax 



approaches for t 



the function 



-7 fl - 52N(e-3)/2 

1™ P{^, t)^- ^ „.Ji 1, 0{l - x') , 



(27) 



(28) 



In this case wc have (|24|) as stable fixed point. For (Qo) = 
the average coordinate converges to again, but the 
second moment approaches the e-independent maximum 
value 



lim (X2(t)) 



r 



(29) 



As the probability density for the coordinate P{x, t) has 
support limited to the interval [—g/7, (7/7], it follows from 
(|29|l the the limiting probability density is composed of 
two 5-functions of equal weight ^ located at the edges 
of the latter interval. More generally, for non-symmetric 
initial condition for the noise, (Qo) ^ 0, the limiting prob- 
ability density is the sum of two (5-functions, 



lim P{x, t) = p^S 



q 

X 

7 



P- 



x+^- 

7 



(30) 



the weights of which depend non-trivially on e and the 
initial condition 



I 



P± = 



f ± 



(Qo) 



aie) 



(31) 



where cr(e) = limt^tx) si(t). 

The function A(t) relaxes to zero. Therefore, in this 
regime, the switching of the external force asymptotically 
stops and the coordinate X{t) approaches either the value 
+q/j or —(7/7, where it freezes. So, the coordinate ac- 
quires a random but time-independent asymptotic value. 
More precisely, the mean coordinate (X(t)) approaches a 
generally non-zero asymptotic value, which depends on 
the initial condition. This is the manifestation of glassy 
state in the regime e < 0, characterised by broken er- 
godicity and non-zero Edwards-Anderson order parame- 
ter. This point will be discussed more in detail later in 
the presentation of correlation functions. 

As in the ergodic phase, all quantities relax toward 
their limit values exponentially for large times. The rate 
of convergence is governed by the eigenvalue with smallest 
modulus, which is now 



Ml 



-27(1 -e) 



(32) 



Note that, contrary to the ergodic regime, the eigenvalue 
is always a real number, so no oscillations occur, at least 
in the linearised approximation. 
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3.5 Marginal case 



1 



Let us proceed by approaching the marginal case e = 1 
from the non-crgodic side, i. e. from below. It might be in- 
structive to cast the equations (|22|l in terms of the eigen- 
modes of the linearised approximation. Namely, we can 
introduce the functions 



m 



1 



^^ I -S3(i) + 2e7S4(i) + ^ 
(2e- 1)7 V 7 

and express the equations H22() in the form 

e--27(l-e)e + 2e73(2ee + ??)? 
77 —777 — 2e7'''(2e^ + 77) ^ 



(33) 



(34) 



The function ^(t) has a straightforward physical interpre- 
tation: it describes the time evolution of the switching rate 
of the external force. Indeed, inserting into we 
get 

A(i) = 27^(1 -2e)e(t). (35) 

The equations (|34|l are a convenient starting point for 
the investigation of the marginal regime. Taking e = 1, the 
linear term in the equation for ^ vanishes, while the linear 
term in the equation for 77 remains. This suggests that 
in the long-time regime the value of 77 will be negligible 
compared to ^. This consideration will yield the leading 
term in the relaxation. 

Thus, supposing I77I <C \£,\ we get the following approx- 
imate equation for ^ 



which leads to the following asymptotic behaviour 

1 1 
4"^ 1 



(36) 



(37) 



Now we must check the assumption that 77 is negligible 
compared to ^. However, from H34|) we can see that the 
leading term in the relaxation of 77 is 



/ N 11 



00 



(38) 



and the assumption is therefore consistent. 

The consequence to draw is that in the marginal regime 
the relaxation becomes power- law with exponent —1. Es- 
pecially, the relaxation of the switching rate follows the 
behaviour 

m c^Yt ' * ^ «^ ■ (39) 

In Fig.2|we can compare the numerical solution with the 
asymptotic behaviour (|39|l . We can see not only that that 
the function A(t) approaches zero according to the power 
decay , but also the corrections to the asymptotic be- 
haviour can be well approximated by a power. Indeed, 
from the inset in Fig. 0] we can see that 



2tX(t) 



- 1 ~ 3t 



-0.9 



(40) 
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Fig. 4. Time evolution of the switching rate in the marginal 
regime e = 1, for 7 = 1, q = 1. The red line is the numerical 
solution, blue line the asymptotic analytical solution 1391 . In 
the inset the deviation from the expression (1391 is shown (red 
line). The blue line is power dependence 3 x t"'^'^ . 



It is interesting to note that the power in the correction 
is not an integer, so the naive expansion of the solution in 
powers of cannot be used here. Instead, the behaviour 
()40|l suggests the expression in the form of a continued 
fraction 



(41) 



where the values ai = 2 and ai = 1 are known exactly 
and the next pair of parameters is estimated from the 
numerical solution as a2 ~ 1/6 and 02 — 0.9. 



4 Correlation functions 

Additional information on the properties of the transition 
from ergodic to non-ergodic behaviour which occurs at the 
value e = 1 can be gained from the two-time correlation 
functions. Let us have t > ti > and define the correlation 
function 

C(t,<i) = (X(i)X(ti)) . (42) 

It can be expressed through the functions si(i) to 34(1). 
The most general formula is 

C(t,ii) = (Xo)e-^(*+*i) + 

+ (XoQo) [e-^*^ S2it)+e-^'s2ih)] + 

+ 2e-'^(*-*i) S4(ti)+ 

53(^1 



(43) 



S2(t)-e-''(*-*i)s2(ii 



although we suppose throughout this section that Xq = 0. 
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(X(t)X(ti)) 




Fig. 5. Correlation function in the ergodic regime, for e = 1.2, 
7 = 1, g = 1. 



We show in figures (ergodic regime) and El (non- 
ergodic regime) the evolution of correlation functions us- 
ing the numerical solution for Si to S4. Wc can observe the 
damping of the correlations in the ergodic regime, while 
in the non-ergodic regime the correlations converge to a 
finite limit. Let us now turn to the analytic investigation 
of the long-time behaviour of the correlation function. 

For long enough times we can suppose that wc are 
in the regime of exponential asymptotic relaxation of the 
functions si{t) to S4{t) and X(t), which is governed by the 
eigenvalue closest to 0, as given by (|26|) and (|32|l . 

Let us start with the ergodic regime e > 1. We find 
that for both i — > 00 and ti s- 00 the correlation function 
behaves like 



7M2-e) 



-7(e-l)(t-ti)_ 



-(e-l)e 



-7(t-*l) 



(44) 



Important feature of this result is that the correlation in 
asymptotic regime depends only on the time difference 
t — ti and decays to zero when this difference increases. 
This supports the picture of the e > 1 phase as a usual 
ergodic regime without long-time correlations. 

The situation is dramatically different for e < 1. Tech- 
nically speaking, it is important that the function A{t) has 
a finite limit for large times, limt^oo A{t) = ^00 < 00. For 
e close to 1 (i. e. 1 — e ^ 1) the approach to this limit 
value can be written in the form 



(45) 



where is a constant depending on e. 

For large times the correction will be small and we 
can formally write the correlation function as expansion 
in powers of 9. However, we should bear in mind that it is 
not d itself, which is small, but the factor e^^^^^^*^)* which 
appears always together with 6. 




Fig. 6. Correlation function in the non-ergodic regime, for 
6 = 0.8, 7 = 1, q=L 



Finally, 



72 (7 - 
(7 - ^)(27 - 



(46) 



using /i = 27(1 — e) for shorter notation. 

We can see that the qualitative difference from the 
ergodic regime consists in the fact that for e < 1 the 
correlation function converges to a positive e-independent 
constant q^/7^. If we define the Edwards- Anderson order 
parameter 



lim lim 

T — >00 tl — >0 



C{h+T,tl) 



(47) 



we can see that qea jumps discontinuously from the value 
qea = for e > 1 to qea ~ q'^/j^ for e < 1. This observa- 
tion represents another evidence that there is a transition 
from ergodic regime to non-ergodic glassy regime at e = 1 . 
As the Edwards- Anderson parameter is discontinuous at 
the critical point, the transition should be classified as 
first-order from this point of view. However, because we 
do not deal with an equilibrium transition and the phe- 
nomenon is of purely dynamical origin, the canonic clas- 
sification of phase transition as first or second order is of 
limited relevance here. 



5 Critical behaviour at e ^ 1 

We have already seen that it is possible to characterise the 
glass transition at e = 1 through the Edwards-Anderson 
order parameter qea- It has discontinuity at the transi- 
tion, so the corresponding critical exponent is 0. Here we 
investigate another quantity, which can play the role of an 
order parameter, being zero in the ergodic and non-zero 
in the non-ergodic phase. 
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Fig. 7. The critical behavior at e — > 1. Points (+) are results 
of numerical integration. The green lino is the approximation 
l|55|l . while the blue line is the exact upper bound crupper (e) 
given by formula 11541 . 

The quantity in question will describe how the initial 
conditions affect the asymptotic vale of the average coordi- 
nate. We have already touched this point in Sec. 13. 41 Up to 
now we assumed that the initial condition for the noise is 
such that (Qo) = 0. This assumption has the consequence 
that both in ergodic and non-crgodic regime the average 
coordinate converges to 0. In this section we investigate 
the case (Qo) ^ 0. 

From (|15|1 and (|18|) we can see that 

lim(X(<)) = Ma(6) (48) 

where we defined, as in Sec. 13.41 

cr(e) ^ lim si{t) (49) 

t — >oo 

stressing explicitly the dependence on e. The equations 
hold also in case (Qo) 7^ 0. Therefore, we can proceed 
without further complications. We only need to plug the 
solution obtained the same way as in the sections 13 and 0] 
into the definition I|1HI of the function si(t) and find the 
Umit igni- 

Let us first present the results of numerical integration. 
We will turn to the analytic estimate afterwards. The Fig. 
13 shows the results of numerical integration, indicating 
that asymptotically for e ^ 1 the behaviour follows the 
power law 

aie)^il-e)i. (50) 



set of equations for the pair tp{t) and X{t) 

''7 (51) 
A = -27(l-e)A + 2e72 XiP 

with initial conditions ijj{0) = — 1/7, A(0) = ej. If we 
further define Aoo = /g X{t)dt and iFoo = Jq '4'{t)dt, we 
can integrate both LHS and RHS of the equations (|51|l 
and obtain the exact relation between Aoo and if'oo 

£ = -2Aoo-2e-f^^oo ■ (52) 

Knowing Aoo would solve the problem, because (T(e) = 
e~"^°°. However, in addition to (|52|) we need some other 
condition. It can be established from the observation that 
the equation for X{t) can be formally solved in the form 

A(t) = £7 cxp ^-27(1 - e)t + 2e72^ V(i')dt'^ • (53) 

and because ip{t) < 0, we have J^ip{t')dt' > ^oo- There- 
fore, we can write the following upper bound 

a{e) < exp (^-^Wl (f^)) = ^"pp-(^) (^4) 

where Wl {x) is the Lambert function defined by the equa- 
tion Wi^{x) e'^i'(^) = X. 

The leading term in the asymptotic behaviour of the 
Lambert function for large argument is Wi^{x) ~ Incc. If 
we use it as an approximation for calculating the asymp- 
totic behaviour of a{e), starting with H54|) we finally get 

CT(e)~Ve\/T^ (55) 

which is compatible with the behaviour (|50|l . Actually, we 
can see in the Fig. [3 that the approximation H55(l fits very 
well the results from numerical integration and lies much 
closer than the exact upper bound H54() . Thus, we conjec- 
ture that the formula l|55(l is in fact the correct asymptotic 
behaviour for e — > 1. 

To sum up, in the regime e > 1, the initial conditions 
are irrelevant for long-time dynamics, as expected in the 
ergodic phase. On the other hand, for e < 1, we observe 
that the asymptotic value of the average coordinate de- 
pends on the initial condition for the noise, which is yet 
another signature of ergodicity breaking. The factor a{e) 
measures the sensitivity to initial conditions: it vanishes in 
ergodic phase but remains non-zero in non-crgodic phase. 
So, it may be considered as a kind of order parameter. 
Close to the critical point e = 1 it approaches zero contin- 
uously as a power with critical exponent 1/2. 



The exponent 1/2 is still observed only empirically and 

we do not possess any proof that this is the exact value. 6 Response tO harmonic perturbation 

However, we may get some analytical argument in favour 
of this type of behaviour from the equations (|19|l . (|20(l . Let us investigate now the response of the particle to 
and H21|) . which can be rewritten in slightlv different form, the external driving force. Adding the additional term 
Defining new function i/j{t) = S3(i) — I/7 we can write the F{t) = Fq cos{ujt) at the right hand side of Eq. one 
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Fig. 8. Response of the switching rate to harmonic external 
driving force for 7 = 1, g = 1, e = 0.8, uj = 2-n and Fq = 0.7 
(blue line) and Fq — 0.15 (green line). The red line shows the 
time dependence in absence of the external driving force. 




10-^ 10-^ 10-1 1 



i^o 

Fig. 9. Dependence of the constant term in the Fourier series 
15911 for \st{t) on the amplitude of the driving force (+), in the 
regime e < 1. The parameters are 7 = 1, g = 1, e = 0.8 and 
the frequency is = 7r/2. Full line is the function 1.4 Fq^- 



can repeat, mutatis mutandis, all steps leading to the sys- 
tem of equations for the functions si{t) to S4{t). We find 
that the equations lfT7|) - (pn|) hold unchanged, while the in- 
fluence of the external force modifies the expression H21|l 
for X{t). Actually, in the present case one gets 



m 

£7 



1-2753(0+272 Si{t) + 

^ (^-le-^' j\{t')e"' dt'-F{t) 



(56) 



We assumed (Qo) = and (Xq) = here. 

First quantity to study is the response of the average 
coordinate. We find 



e"^* / F(t') C^*' dt' 
"'0 



(57) 



Obviously enough, in the stationary regime the average 
coordinate oscillates around with the same frequency as 
the driving force F{t). We arrive at the standard Debye- 
type dynamic susceptibility 



x{t - 1') = eit - t')e-^^'-*'^ 



(58) 



Thus, the exact response in terms of the coordinate is 
linear. This behaviour also does not depend on the value 
of e. 

The situation becomes much more complicated when 
we turn to quantities, which depend non-linearly on the 
coordinate, especially the switching rate X{t). We solved 
numerically the set of equations H17() to (|20|l and (|56(l . 
We can see in Fig. |S1 the evolution of the function A(t) 
within the non-crgodic regime, with e = 0.8. Comparing 



the behaviour with the evolution in absence of the exter- 
nal harmonic perturbation we can see that the switching 
rate X{t) does not approach zero any more, but it oscillates 
around some finite value, which we will denote ^o- This 
qualitative feature holds for whatever small external force. 
However, when the amplitude Fq of the external field goes 
to zero, also the value of Aq approaches zero according to 
Aq ~ Fq^, as can be seen from Fig.^ Generalising the lin- 
ear stability analysis of Sec. 18.21 to harmonic oscillations, 
we conclude that the external field continuously shifts the 
fixed point with A = 0, (i. e. also Aq = 0) to a position 
with positive Aq, but the value of the shift vanishes when 
the amplitude of the perturbation goes to zero. 

Figure ^| exemplifies the evolution of the functions 
33(1) and S4{t). We can clearly see that the oscillations are 
not harmonic. Generally, in the stationary regime these 
functions are non-harmonic but periodic with the doubled 
frequency 2a;. Thus, the same holds also for the function 
X{t). Using H56|) the stationary response in terms of the 
switching rate can be written as Fourier series 



t (t) = + ^ (Ak sin 2kujt + Bk cos 2kujt) 



(59) 



The amplitudes of the harmonic components Aq, Ak,Bk, 
k = 1,2, satisfy a complicated infinite set of quadratic 
equations. 

To asses the weight of the higher harmonics we per- 
formed the fast Fourier transform of the time evolutions 
obtained by numerical solution, throwing away the initial 
transient regime. To illustrate the presence of higher har- 
monics we chose the function S3(<). The modulus of its 
Fourier transform 83(1') = J 53(^)6"^'^"^* dt is shown in 
Fig. ^3 We can clearly see the peaks at the multiples of the 
basic frequency. We can also observe that the higher har- 
monics have quite considerable weight. In the inset of Fig. 
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Fig. 10. Response of the functions S3{t) (red line) and S4{t) 
(blue line) to harmonic external driving force for 7 = 1, q = 1, 
e = 0.8, = TT and Fo = 10. 



Illlwe show also the Fourier transform of the function X{t). 
Here, the higher harmonics are much less pronounced. 

The most important feature of the time evolution un- 
der the influence of external harmonic force is the obser- 
vation that the switching rate remains always positive. 
This leads to the already mentioned fact that whatever is 
the coupling strength parameter e, the mean coordinate 
in stationary regime oscillates around zero, irrespectively 
of the initial conditions. However, this is the signature of 
ergodicity, so the glassy behaviour disappears under the 
influence of arbitrarily small external perturbation. Such 
a behaviour was already observed also in the model of 
sheared colloid \2B\ . 

Both observations can be easily understood. The ex- 
ternal force drags the particle back and forth. Once mov- 
ing, the particle induces through the back-reaction the 
fluctuations of the environmental force, which prevents the 
system from freezing in a non-crgodic state. This holds for 
any positive amplitude of the driving force. However, when 
the force diminishes, there is still longer transient period, 
where the system apparently relaxes toward the arrested 
state, as can be seen qualitatively in Fig.|Hl In the limit of 
infinitesimally small driving, the transient time blows up 
and the asymptotic state corresponds to the dynamically 
arrested state. This picture is consistent with the view of 
glassiness purely dynamical phenomenon. 

We also observed the response of the system to a signal, 
which is switched on only after the system relaxed very 
close to the arrested state. The results can be seen in Fig. 
1121 Initially, the switching rate relaxes toward zero, but 
after the perturbation it settles on oscillating behaviour. 
The average coordinate initially approaches non-zero value 
(we have chosen initial condition (Qo) > 0), but the per- 
turbation brings it to oscillations around zero. This can 
be interpreted as a schematic picture of shear thinning, 
although the model is too much simplified to account for 
the shear thinning quantitatively. 



Fig. 11. Fourier transform of the stationary oscillations of the 
function 53(4) under the influence of harmonic external driving 
force for e = 0.8, uj = n and Fo = 10. It was calculated using 
the fast Fourier transform algorithm. In the inset the Fourier 
transform of oscillations of the switching rate. The finite width 
of the peaks, the noise and the continous part of the spectrum 
in the inset are due to numerical imprecision of the fast Fourier 
transform procedure. 



7 Conclusions 

Motivated by the dynamical arrest phenomenon in colloids 
we formulated and solved a stochastic dynamical model, 
where the coordinate of a single particle evolves under 
the influence of stochastic environmental force. The back- 
reaction couples the switching rate of the force to the av- 
erage of the square of the velocity of the particle. The 
strength of the coupling e is the crucial parameter which 
determines the behaviour of the system. The problem re- 
duces to a set of coupled non-linear differential equations, 
which was investigated both analytically and numerically. 

The back-reaction induces a phase transition from the 
ergodic phase for e > 1 to the non-ergodic glassy state 
for e < 1. The transition is observed qualitatively in the 
behaviour of the switching rate, decaying to zero in non- 
ergodic state, while staying positive in the ergodic state. 
The Edwards- Anderson parameter, established from the 
two-time correlation functions, is discontinuous at the tran- 
sition; it is zero in ergodic phase (e > 1), while for e < 1 
it acquires finite value independent of e. The critical point 
e = 1 is characterised by power-law decay of the switching 
rate. The leading term ^ in the long-time behaviour 
was calculated analytically. 

We investigated the critical behaviour at e ^ 1 through 
the dependence of the average coordinate in the long-time 
limit on the initial condition. We find that the asymptotic 
value of the average coordinate is proportional to the aver- 
age initial value of the force, where the proportionality fac- 
tor (T(e) is singular at the transition. In the ergodic phase 
we have a{e) — identically, while in the non-crgodic 
phase, close to the transition, we found cr(e) ~ (1 — e)^/^ 
for e ^ 1~. 
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Fig. 12. Time evolution of the verage coordinate (a) and 
switching rate (b) in the regime e < 1, when the external 
harmonic perturbation was switched on at time t = 50. The 
parameters are ^ = 1, q = 1, e — 0.8, the amplitude of the ex- 
ternal perturbation is Fo = 0.3 and the frequency is lj = tt/2. 
The initial condition is (Qo) = 1- 



Therefore, we find a situation quite unusual from the 
point of view of static equilibrium phase transition. In- 
deed, we have two variables, which may be considered as 
order parameter, namely the Edwards- Anderson parame- 
ter and the quantity cr(e). While the former is discontin- 
uous at the critical point, thus indicating first-order tran- 
sition, the latter is continuous, suggesting second-order 
transition, the discrepancy is to be attributed to purely 
dynamical nature of the transition. 

Finally, we investigated the response of the system to 
harmonic external perturbation. We found that the exact 
response of the coordinate is linear. On the other hand, 
the response of the variables which are quadratic functions 
of the coordinate and velocity, like the switching rate, was 
non-linear and generically contains all higher harmonics, 
as was seen in the Fourier transform of the signal. We 
also observed that arbitrarily weak external perturbation 
is sufficient to "melt" the non-ergodic glassy state and 
bring it back to ergodic behaviour. We may relate this 



feature to the notion of stochastic stability "TT; in this 
view our system is not stochastically stable. However, our 
finding is in accord with previously observed behaviour of 
sheared colloids [23 ■ It makes also connection to the rhe- 
ological properties of thixotropic fluids |48ll49j . although 
our model is too simplified to give quantitative predictions 
in this direction. 

The back-reaction mechanism described by (|5J repre- 
sent the simplest choice. One may ask what would hap- 
pen if we tried another prescription. We expect that the 
methods used here will be as well applicable if we gen- 
eralise (0) as X{t) = (i^(V^(t))) for an analytic function 
F{x). More complicated situation would appear if the de- 
pendence was non-local in time, e. g. of the form X{t) = 
/*(V(t)V(t'))) K{t-t') At' with some kernel K{t). Such an 
approach would bring our model closer to the well-studied 
Mode Coupling equations, but it goes beyond the scope 
of the present work. 

We want to dedicate this paper to the memory of Prof. Vladislav 
Capek, a passionate theoretical physicist, our dedicated teacher 
and good friend, who passed away shortly before this work was 
completed. 

This work was supported by the project No. 202/00/1187 
of the Grant Agency of the Czech Republic. 



References 

1. C. A. Angell, Science 267, 1924 (1995). 

2. R. G. Palmer, D. L. Stein, E. Abrahams and P. W. An- 
derson, Phys. Rev. Lett. 53, 958 (1984). 

3. G. Paladin, M. Mezard and C. de Dominicis, J. Physique 
Lett. 46, L-985 (1985). 

4. T. R. Kirkpatrick and D. Thirumalai, Phys. Rev. Lett. 58, 
2091 (1987). 

5. T. R. Kirkpatrick and D. Thirumalai, Phys. Rev. B 36, 
5388 (1987). 

6. T. R. Kirkpatrick, D. Thirumalai and P.G. Wolynes, Phys. 
Rev. A 40, 1045 (1989). 

7. F. Ritort, Phys. Rev. Lett 75, 1190 (1995). 

8. S. Franz and J. Hertz, Phys. Rev. Lett. 74, 2114 (1995). 

9. L. F. Cugliandolo and J. Kurchan, Phys. Rev. Lett. 71, 
173 (1993). 

10. S. Franz, M. Mezard, G. Parisi, and L. Peliti, Phys. Rev. 
Lett. 81, 1758 (1998). 

11. A. Barrat, R. Burioni and M. Mezard, J. Phys. A: Math. 
Gen. 29, 1311 (1996). 

12. J. Kurchan, J. Phys. I France 2, 1333 (1992). 

13. L.F. Cughandolo and J. Kurchan, Phil. Mag. B 71, 501 
(1995). 

14. E. Marinari, G. Parisi and F. Ritort, J. Phys. A: Math. 
Gen. 27, 7647 (1994). 

15. T. M. Nieuwenhuizen, J. Phys. A: Math. Gen. 31, L201 
(1998). 

16. A. Cavagna, I. Giardina, and G. Parisi, Phys. Rev. Lett. 
83, 108 (1999). 

17. A. Cavagna, I. Giardina, and G. Parisi, J. Phys.: Condens. 
Matter 12, 6295 (2000). 

18. M. Mezard and G. Parisi, Phys. Rev. Lett. 82, 747 (1999). 



12 



Frantisek Slanina, Petr Chvosta: Glass transition in a simple stochasic model with back-reaction 



19. B. Coluzzi, G. Parisi, and P. Verocchio, Phys. Rev. Lett. 
84, 306 (2000). 

20. L.F. Cugliandolo, J. Kurchan, G. Parisi and F. Ritort, 
Phys. Rev. Lett. 74, 1012 (1995). 

21. G. Parisi and F. Slanina, Phys. Rev. E 62, 6554 (2000). 

22. K. Kawasaki and B. Kim, Phys. Rev. Lett. 86, 3582 

(2001) . 

23. S. M. Fielding, P. SoUich, and M. E. Gates, J. Rheol. 44, 
323 (2000). 

24. M. Fuchs and M. E. Gates, Faraday Discussions 123, 267 
(2003). 

25. M. E. Gates, preprint 'cond-mat /02 1 10661 

26. A. M. Puertas, M. Fuchs, and M. E. Gates, Phys. Rev. 
Lett. 88, 098301 (2002). 

27. M. Fuchs and M. E. Gates, Phys. Rev. Lett. 89, 248304 

(2002) . 

28. E. Zaccarelh, G. Foffi, K. A. Dawson, S. V. Buldyrev, 

F. Sciortino, and P. Tartagha, Phys. Rev. E 66, 041402 
(2002). 

29. A. Lawlor, D. Reagan, G. D. McGuUagh, P. De Gregorio,l 
P. Tartaglia, and K. A. Dawson, Phys. Rev. Lett. 89, 
245503 (2002). 

30. E. Stiakakis, D. Vlassopoulos, B. Loppinet, J. Roovers, and 

G. Meier, Phys. Rev. E 66 , 051804 (2002). 

31. L. Fabbian, W. Gtze, F. Sciortino, P. Tartaglia, and F. 
Thiery , Phys. Rev. E 59, R1347 (1999). 

32. M. E. Gates, J. P. Wittmer, J.-P. Bouchaud, and P. 
Glaudin, Physica A 263, 354 (1999). 

33. E. Bertrand, J. Bibette, and V. Schmitt, Phys. Rev. E 66, 
060401 (2002). 

34. W. Gotze, in: Liquids, Freezing and Glass Transition, ed. 
J. P. Hansen, D. Levesque and J. Zinn- Justin, (North- 
Holland, Amsterdam, 1991) p. 287. 

35. W. Gotze and L. Sjogren, Rep. Prog. Phys. 55, 241 (1992). 

36. J. P. Bouchaud, L. F. Gugliandolo, J. Kurchan, M. and 
Mezard, in: Spin glasses and random fields, ed. A. P. 
Young, (World Scientific, Singapore, 1997) p. 161. 

37. L. Sjogren, Physica A 322, 81 (2003). 

38. L. Gughandolo and J. Kurchan, Physica A 263, 242 (1999). 

39. P. Chvosta and F. Slanina, J. Phys. A: Math. Gen. 35, 
L277 (2002). 

40. W. Horsthemke and R. Lefever Noise-Induced Transitions. 
Theory and Applications in Physics, Chemistry and Bi- 
ology Springer Series in Synergetics, vol. 15, (Springer, 
Berlin, 1984). 

41. D. R. Cox Renewal Theory (Wiley, New York, 1962). 

42. D. L. Snyder Random Point Processes (Wiley, New York, 
1975). 

43. Y. Pomeau J. Stat. Phys. 24, 189 (1981). 

44. P. Sibani and N. G. van Kampen Physica A 122, 397 
(1983). 

45. A. Morita Phys. Rev. A 41, 754 (1990). 

46. M. Abramowitz and 1. A. Stegun (editors) Handbook of 
Mathematical Functions (Dover, New York, 1970). 

47. G. Parisi, preprint [cond-mat/0007347| 

48. G. Voltz, M. Nitschke, L. Heymann, and I. Rehberg, Phys. 
Rev. E 65, 051402 (2002). 

49. P. Coussot, Q. D. Nguyen, H. T. Huynh, and D. Bonn, 
Phys. Rev. Lett. 88, 175501 (2002). 



